function [etrh, detrb, etri, etr2i] = hydro(filename, metricorder)

  load(filename);


  % summary
%   disp(sprintf('Grid Size   : (%d, %d, %d)', params.Nx, params.Ny, params.Nz));
%   disp(sprintf('Grid Length : %f', params.Lx));
%   disp(sprintf('Cell Length : %f', params.Dx));
%   disp(sprintf('Target Time : %f', params.Tfinal));
%   disp(sprintf('Time steps  : %d',  params.Tsteps));
%   disp(sprintf('Hydro Size  : %6.2fMB', 3*params.Nx.*params.Ny.*5.*TIME_STEPS.*64/8/1024/1024));
%   disp(sprintf('Metric Size : %6.2fMB', length(params.metricEvaluationPoints).*params.Nx.*params.Ny.*params.Nz.*4.*4.*3.*64./8./1024./1024.*(params.timeOrder+1)));

  % prepare figures for plotting
  close all;

  % load color maps
  load('cm', 'mycmap');
  params.colormapOffset = 512;

  if (params.plotHydro || (params.plotMetric && params.metricOn))
    fig  = figure; 
    set(fig, 'Colormap', mycmap);
    set(fig, 'WindowStyle', 'docked');
  end
  refresh();
  pause(0.001);


  % initialize derivative kernels
  h2dinit(params);
  evalpoints = params.metricEvaluationPoints;

  etrh   = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  detrb  = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  etri   = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  etr2i  = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  

  for i=1:length(metricorder)
    params.metricOrder = metricorder(i);
    [g gt gtt] = h2metric(params, tOut, uOut, udtOut, udtdtOut);
    tic
    [etrh(i,:,:,:), detrb(i,:,:,:), etri(i,:,:,:), etr2i(i,:,:,:)] = h2gr(params, g, gt, gtt);
    disp(sprintf('h2gr in %.2f sec', toc))
  end

end


















